
#include "Data.h"
#include "ICs.h"
#include "Kernel.h"

using namespace std;

class MyICs : public ICs {
	public:
		
		MyICs(Data &data): ICs(data) {}
		
		int nfloors;
		
		void make_WTC(){
	
			grid.xn = 2*64;
			grid.yn = 2*64;
			grid.zn = 300;
	
			grid.xs = 3;
			grid.ys = 3;
			grid.zs = 3;
	
			data.pt.resize(2);
	
			data.pt[0].m = 1;
			data.pt[1].m = 1;
	
			data.pw[0][0].s  = 1;
			data.pw[0][0].rc1f = 1;
			data.pw[0][0].rc2 = 2;//pow(1.5, 2);//SQR(2);
	
	
			data.pw[1][1].s  = 1.75;
			data.pw[1][1].rc1f = 1;
			data.pw[1][1].rc2 = 2.5;//SQR(2);
	
			data.pw[0][1].s  = 1;
			data.pw[0][1].rc1f = 1;
			data.pw[0][1].rc2 = 1.7;//pow(1.2, 2);//SQR(2);	

			data.mU = 1; // kg
			data.lU = 1.0/6; // meters
			data.tU = 1; //sec
			data.vU = data.lU/data.tU;
	
			double lU = data.lU;
			double vU = data.vU;
			
			// core columns
			
			addHCPBox(0, 1, dVector( -15,0,4*32-25), Vector<int>(32,32,32), dVector(1,1,1), dVector(0,0,-100), 0);
			addHCPBox(0, 1, dVector( -15,0,3*32-25), Vector<int>(32,32,32), dVector(1,1,1), dVector(0,0,-100), 0);
			addHCPBox(0, 1, dVector( -15,0,2*32-25), Vector<int>(32,32,32), dVector(1,1,1), dVector(0,0,-100), 0);
			addHCPBox(0, 1, dVector( -15,0,1*32-25), Vector<int>(32,32,32), dVector(1,1,1), dVector(0,0,-100), 0);
			
			
			remove_overlaps();
		
			data.pw[0][0].e  = (500000.0/12)*SQR(data.tU)/SQR(data.lU); 
			data.pw[1][1].e  =      (1.0/12)*SQR(data.tU)/SQR(data.lU);
			data.pw[0][1].e  = (500000.0/12)*SQR(data.tU)/SQR(data.lU);

			data.update_derivatives();
			
		}
	
};
	


int main(int argc, char *argv[]) {

	
	Data data;
	data.tU = 0.001; // s
	
	
	MyICs ics(data);
	
	ics.kfac = 0.985;
	ics.nfloors = 15;  
	ics.make_WTC();
	ics.print();
	
	
	
	Kernel kernel(data);
	
	kernel.period = 10; // s
	kernel.ss_period = 0.01; // s
	
	kernel.grid.xs = 3;
	kernel.grid.ys = 3;
	kernel.grid.zs = 3;
	
	kernel.grid.xn = 2*75;
	kernel.grid.yn = 2*75;
	kernel.grid.zn = 300;
	
	kernel.dtmax = 4e-5/data.tU;
	kernel.dRmax = 1.0/10;
	kernel.dEmax = data.pw[0][1].e12/5; // minimum potential well depth 
	
	kernel.outputfolder = "out.8.2";
	
	kernel.relaxationF = true;
	kernel.Frx = 1e-4;
	kernel.Vcl = 0.991;

	kernel.cdt = 1.1;
	kernel.run();
	

	return 0;
	
}
